During the summer of 2021, the Charlottesville Office of Sustainability organized volunteers to collect temperature data from bikes and automobiles on August 24th. These data are now available for download on the City’s website and from the OSF repository. The details of the data collection are available in this report and the machine learning methods used to generatethe rasters are described in these publications:

In addition to the interpolated rasters produced by CAPA Strategies in partnership with NOAA, we can access Landsat images via the USGS EarthExplorer repository and use the thermal bands to estimate land surface temperatures.

A Landsat 8 image was acquired at 3:53 PM on Tuesday August 24, 2021 with low cloud cover (0.31 percent) and cloud shadow (0.44 percent). This is interesting because it coincides precisely with the Heat Watch data collection effort.


As shown in the screen capture above, this image covers our entire study area because the rectangle is the bounding box of the Meadow Creek watershed boundary plus a distance buffer of 1 mile.

library(sf)
library(raster)
library(tidyverse)
library(ggpubr)
library(RColorBrewer)


pm_rast_0 <- raster::raster("./data/rasters_chw_charlottesville_010622/pm_t_f.tif")
landsat_rast_0 <- raster::raster("./data/LC08_CU_027010_20210824_20210904_02_ST_B10.tif")
landsat_rast_1 <- raster::projectRaster(landsat_rast_0, crs=crs(pm_rast_0))
watershed <- st_read("./data/meadow_creek_watershed_boundary.shp")
## Reading layer `meadow_creek_watershed_boundary' from data source 
##   `C:\Now\Collaborations\3Cavs\Data\Heat_Exposure\data\meadow_creek_watershed_boundary.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 11476070 ymin: 3899273 xmax: 11492140 ymax: 3919360
## Projected CRS: NAD83 / Virginia South (ftUS)
bounding_box_0 <- st_read("./data/bbox_meters.shp")
## Reading layer `bbox_meters' from data source 
##   `C:\Now\Collaborations\3Cavs\Data\Heat_Exposure\data\bbox_meters.shp' 
##   using driver `ESRI Shapefile'
## replacing null geometries with empty geometries
## Simple feature collection with 1 feature and 10 fields (with 1 geometry empty)
## Geometry type: GEOMETRYCOLLECTION
## Dimension:     XY
## Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
## Projected CRS: WGS 84 / UTM zone 17N
bounding_box <- st_transform(bounding_box_0, st_crs(pm_rast_0))

pm_rast_spdf <- as(pm_rast_0, "SpatialPixelsDataFrame")
pm_rast_df <- as.data.frame(pm_rast_spdf)
colnames(pm_rast_df) <- c("value", "x", "y")
wshed_capa <- st_transform(watershed, st_crs(pm_rast_0))

capa_heat <- ggplot() +
  geom_raster(data = pm_rast_df, aes(x = x, y = y, fill = value)) +
  scale_fill_gradientn(colors = brewer.pal(n = 8, name = "RdYlBu"), trans = "reverse") +
  labs(fill = "Temp. °F", title ="Late Afternoon Temperature - CAPA Heat Watch", subtitle = "24 August 2021") +
  geom_sf(data = wshed_capa, colour = "lightgreen", fill = NA) + 
  theme_void()
capa_heat

# Crop to reasonable spatial extent (not entire scene)
wshed_lst <- st_transform(watershed, st_crs(landsat_rast_1))
bbox_extent <- extent(711493, 726111, 4209536, 4220981)
landsat_rast_2 <- raster::crop(landsat_rast_1, bbox_extent)


# Convert to Fahrenheit after scaling raw values to Kelvin
# per this reference   https://www.usgs.gov/faqs/how-do-i-use-scale-factor-landsat-level-2-science-products
landsat_rast_K <- ((landsat_rast_2 *  0.00341802) + 149.0)
landsat_rast_F <- (((landsat_rast_K - 273.15) * 1.8) + 32)


landsat_rast_3 <- crop(landsat_rast_F, extent(pm_rast_0)) 

# Convert to an object ggplot2 can handle
landsat_rast_spdf <- as(landsat_rast_3, "SpatialPixelsDataFrame")
landsat_rast_df <- as.data.frame(landsat_rast_spdf)
colnames(landsat_rast_df) <- c("value", "x", "y")


landsat_lst <- ggplot() + 
  geom_raster(data = landsat_rast_df, aes(x = x, y = y, fill = value)) +
#  scale_fill_viridis_c(option = "plasma") + 
  scale_fill_gradientn(colors = brewer.pal(n = 8, name = "RdYlBu"), trans = "reverse") +
  labs(fill = "Temp. °F", title ="Land Surface Temperature", subtitle = "24 August 2021 @ 3:53 PM") +
  geom_sf(data = wshed_lst, colour = "lightgreen", fill = NA) + 
#  coord_sf(ylim = c(layer_scales(capa_heat)$y$range$range), xlim = c(layer_scales(capa_heat)$x$range$range)) +
  theme_void()
landsat_lst

ggarrange(capa_heat, NULL,  landsat_lst, ncol = 3, nrow = 1, widths = c(1, 0.5, 1))

ggsave("capa_temp_and_landsat_lst_24_aug_2021.png", width = 24, height = 16)


You can see above that the CAPA Heat Watch rasters do not cover the entire Meadow Creek watershed.

We can also get the historical weather data from this site collected at the Charlottesville Airport for context.


bbox_extent <- extent(711493, 726111, 4209536, 4220981)
landsat_rast_F_cropped <- crop(landsat_rast_F, bbox_extent)

# Convert to an object ggplot2 can handle
landsat_rast_spdf <- as(landsat_rast_F_cropped, "SpatialPixelsDataFrame")
landsat_rast_df <- as.data.frame(landsat_rast_spdf)
colnames(landsat_rast_df) <- c("value", "x", "y")

ggplot() + 
  geom_raster(data = landsat_rast_df, aes(x = x, y = y, fill = value)) +
  scale_fill_gradientn(colors = brewer.pal(n = 8, name = "RdYlBu"), trans = "reverse") +
  labs(fill = "Temp. °F", title ="Land Surface Temperature - Landsat Thermal Band", subtitle = "24 August 2021 @ 3:53 PM") +
  geom_sf(data = wshed_lst, colour = "lightgreen", fill = NA) + 
  theme_void()


Because the Landsat scenes are provided in a custom map project, we get the funny slant shown above (not an issue though).